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Abstract. We report work towards a relativistic formulation for modeling strongly magnetized 
neutron stars, rotating or in a close circular orbit around another neutron star or black 
hole, under the approximations of helical symmc^try and ideal MHD. The quasi-stationary 
evolution is governed by the first law of thermodynamics for helically symmetric systems, which 
is generalized to include magnetic fields. The formulation involves an iterative scheme for 
solving the Einstein-Maxwell and relativistic MHD-Eulcr equations numerically. The resulting 
configurations for binary systems could be used as self-consistent initial data for studying their 
inspiral and merger. 



1. Introduction 

A uniformly rotating neutron star is significantly deformed when the ratio of kinetic energy T 
to gravitational energy W becomes ~ 0.1. When the magnetic field energy, 



becomes such that A^/|H^| 0.01, the contribution of the magnetic field to the structure of the 
neutron srar may not be neglected. This is estimated to occur for 



Recent observations of anomalous X-ray pulsars, or soft 7-ray repeaters suggest that the neutron 
stars in these systems may be associated with strong mag netic fields around 10^^ - 10^^ G at 
the surface (see, e.g. [1]). Wc may expect that the interior magnetic field of such a strongly 
magnetized neutron star, a magnetar, is a few orders of magnitudes stronger. 

Such strong magnetic fields have not been found in binary neutron star systems, and may not 
survive until merger. Hypothetically, however, strongly magnetized neutron stars or black holes 
may form binary neutron star or black hole - neutron star systems. For instance, the poloidal 








c(t) = V-t o Xt{P) 



Figure 1. Left panel: the history of a sphere S generated by the hehcal vector A;". Right panel: 
a domain of integration S and its boundaries. 

field seen as the surface magnetic field decays but the toroidal field stays strong enough to affect 
the structure of the compact objects. 

In several numerical relativity simulations of magnetized binary neutron stars, it has been 
shown that, after binary inspiral and merger, the magnetic fields are amplified around 10 times by 
magnetic winding and the magnetorotational instability in post-merger, pre-collapse objects [2]. 
Such an object is likely to form a black hole and a magnetized toroid system, which becomes the 
source of a short 7-ray burst. It is desirable to prepare realistic initial data sets for such merger 
simulations calculated by solving the Einstein-Maxwell equations and a first integral of the 
MHD-Euler equation, assuming (quasi-) equilibrium. In this article, we model such magnetized 
binary compact objects in close circular orbits, assuming that the spacetime and magnetic fields 
satisfy helical symmetry and that the stars are in equilibrium [3-5]. 

2. Thermodynamic laws for helically symmetric Einstein-Maxwell spacetimes with 
charged and magnetized perfect fluids 
2.1. Zeroth law 

We consider a globally hyperbolic spacetime {M.,gai3)- The field A;" is transverse to each 
Cauchy surface, not necessarily timelike everywhere, and generates a one-parameter family of 
diffeomorphisms xt- The action of xt on a spacelike sphere 5 on a Cauchy surface generates 
a timelike surface, T{S) = UtXt{<S), called the history of S. Then, k" is a helical vector field 
if there is a smallest T > for which P and Xt{P) are timelike separated for every point P 
outside of the history T{S). A vector /c" written as 

k'' = t'' + Qcf)'' (3) 

is the helical vector, where t"' is a timelike vector and a spacelike vector that has circular 
orbits with a parameter length 27r (see, FUS). In Fig. 1, a schematic figure of the history 7~(5) 
generated by the diffeo xt is presented. 

Each Cauchy surface of a helically symmetric spacetime does not admit flat asymptotics. 
Therefore, the future (past) horizon T-L^ is defined as the boundary of the future(past) domain 
of outer communication of a history T{S) of each spacelike sphere S. If the history T{S) 
of a sphere S is in P^, the future (past) horizon agrees with the chronological past (future) of 
the history T, = dP{T). 

The existence of a global helical symmetry assures that the horizon is a Killing horizon. With 
the null energy condition R^i^l'^l^ > for any null vector the surface gravity k, defined on 
each connected component of the horizon TL^ by 



(4) 



is constant. The proof is given in FUS, in which the conditions of theorems by Priedrich, Racz, 
and Wald [6] are modified to make them suitable for hehcally symmetric spacetimes. Also the 
electric potential in the rotating frame, defined by 

Ea = Fapk^ = -Va$^, = A^k"' + const, (5) 

is constant on Ti^, assuming the symmetry £kAa = 0. Since Ea is null on "H^ and 
Eak'^ = 0, Ea is parallel to the null generator on H^. Hence for any vector 77" tangent to 
H^, r]"'Ea = —rj'^'Va^^ = 0. Therefore, as for the stationary and axisymmetric spacetimes 
shown by Carter [7, 8], the surface gravity k and the electric potential $^ in the rotating frame 
are constant on the horizon H^. 



2.2. First law 

Consider a family of spacetimes, 

Q(A) := [gap{X),u-{X),piX),s{X),Aa{X),r{X)], (6) 
whose Lagrangian density is written 

where n", p, s, e, and are, respectively, the 4-velocity, baryon rest mass density, entropy per 
baryon mass, energy density, and electric 4-current. A generalized first law for the spacetime 
Q(A) associated with the helical symmetry is derived as a variation formula of the Noether 
charge associated with the helical Killing vector defined by [9, 10] 

Q{X) = ^^Q^^dSap, (8) 

"^"^^^ = ^(5"V'-5"V)|a=oV^ 575(A) + ^F^"|a=o[a^(A)-^A^(0)] + 0{X% (9) 

With this choice of 55", the charge Q{X) becomes finite and independent of the sphere S used 
to evaluate the charge Q{X) as long as 5 encloses all black holes and neutron stars. 

To see this, we first check that Q{0) is independent of S. The surface integral of (8) is 
rewritten in terms of integrals over a spacelike hypersurface S transverse to k"" and over the 
black hole boundary Bi, where Bi is the i-th connected component of S H . In Fig.l, the 
region of integration is shown. Since, the boundary of E is the union of the sphere S and Bi, 
dT, = S Ui Bi, a difference Q — J2i Qi calculated at A = 0: 

Q - J] = J^iG'^pk^ - 87rTF^p)dSa " ^ ^ Rk'^dSa- 

+ [■^^■y^'^'^M'" - ^^^AV/3i^"'') dSa-^j^ ^ k^A^F'^f'dSap. (10) 

where Qi{X) is defined on Bi by 

Qi{X) := / Q'^^dSap, (11) 



T^p, the stress-energy tensor of the electromagnetic field, is defined by, 



Tf = ^ (^F-^^F^^ - ^g'^^F^sF''^^ . (12) 

To derive Eq. (10), we have also used a condition 

1 

47r 



^ ^ k^A^F'^^dSa^ = (13) 

on the boundary sphere S to determine the constant of the electric potential in Eq. (5). 
Eq. (10) implies that Q{0) does not depend on the sphere S as long as it encloses all black holes 
and neutron stars, because in the region where the sphere S is located, the integrand of the 
volume integrals in Eq. (10) vanishes when the Einstein and Maxwell equations are satisfied. 

Next, the variation SQ := dQ/d\ in the Noether charge is evaluated in terms of perturbations 
of the baryon mass, entropy, circulation and electric current of each fluid element, and the surface 
areas and charges of the black holes. The detailed calculation is described in [11]: the variation 
dQ becomes 

6Q = f S —^Aispu^'dSa) + ^~^^^ A(pm" dSa) + v'^A{huf)fm''dSa) - Apk^ A{fdSa) 

- irk" - j^k^)AAf, dSa^ + Y, {l^'^i^^ + ^Qf) ' (14) 

where T is the temperature, h is the relativistic specific enthalpy, and f " is the spatial velocity 
defined by the following decomposition of the 4-velocity with respect to the helical vector 

n° = u\k'^ + f°) with v^Vat = and = u'^Vat. (15) 
The electric charge of each black hole Qf is defined by 

Qf:=^<f F-^dSa0, (16) 

which is related to the total electric charge of the system by Stokes' theorem: 

^^'■^hjs ^ L ^ ^ ^^^^ 

Note that the $f , defined on each Bi by 

$f = -A'^ka = + C, (18) 

is constant. In the case of stationary and axisymmetric spacetimes, the mass variation formula 
derived by Carter [7, 8] can be derived from Eq. (14). 

We can now verify that Q{\) is independent of the location of the 2-surfacc S, because the 
charge Q{X) at A = is shown to be independent of 5, and the variation formula (14) implies 
that dQ/dX = SQ is independent of S as long as it encloses the fluid and black holes. Moreover, 
in [11] we have shown that the difference, 

^{Q -T^Qi - ^J^^ k^A^F^^dS^p) =5Q-Y, (^'^i'^A + 5Qf) , (19) 



is invariant under a gauge transformation that respects the symmetry, and hence we verify that 
so is 5Q. 



2.3. Application of the first law to solution sequences in equilibrium 

When inspiraling binary systems, or isolated neutron stars are evolving adiabatically - in a 
timescale much longer than the dynamical timescale - they may be modeled by a sequence of 
solutions in equilibrium. When the first law for a stationary and axisymmetric perfect fluid 
spacetime is applied to a sequence of rotating neutron star solutions whose local changes of 
rest mass and entropy are constant, the first law becomes SM = Q5J. This is a condition for 
applying the turning point theorem [12] to determine the stability of solutions. In this section, 
we consider an application of the first law (14) to a sequence of helically symmetric solutions 
that model inspiraling binary black holes and/or neutron stars, assuming that the neutron star 
matter is a perfect conductor. 

For a perfectly conducting medium, the ideal MHD condition, 

Fa0U^ = 0, (20) 

is satisfied, and hence its curl is written 



i^uFafS = 0, (21) 

which is the magnetic flux conservation law, Alfven's law. For each equilibrium solution, rest 
mass and entropy are conserved: 

£«(PV^) = 0, £us = 0. (22) 



When we consider a sequence of solutions along which the rest mass, entropy, and magnetic flux 
are all conserved, the perturbed conservation laws corresponding to the above Eqs. (21) and 
(22) have first integrals: 

A(pV^) = 0, As = 0, and AF^,^ = 0. (23) 

For black holes, we may assume that areas and charges are constant. With these assumptions, 
the first law (14) is rewritten 

SQ = n v^A{hu^ pu'^ dSa) - Apk^ A{j'^dSa) } . (24) 

Here AF^s = with Fa/s = {dA)ai3, A{dA)ai3 = {dAA)ai3, whence the Poincare lemma implies 
that there exists a function such that AAa = Va*- From a substitution of this to the last 
term in the volume integral of Eq. (14), the latter is shown to vanish. The remaining terms in 
Eq. (24) are related to the circulation of the magnetized flow. 

For a perfect fluid without magnetic fields, the conservation of circulation is written as 
^u^^afs = where UafS •= {d{hu))ai3 is the vorticity tensor. When circulation is conserved along 
a sequence of solutions, it implies that Au^p = 0, and hence the first term in the integral of Eq. 
(24) vanishes [13]. In general, there is no such conservation law for the circulation of magnetized 
flow. However, Bekenstein and Oron [14] (Tarapov and Gorskii [15] for Newtonian MHD) have 
developed a formulation of ideal MHD, in which a generalized circulation of magnetized flow is 
conserved. Their theory is based on a Lagrangian in which the interaction term A^j'^ in Eq. 
(7) is replaced by the term F^fipu^q^ , where the vector q°' is a Lagrange multiplier enforcing 
the perfect conductivity condition F^fju"^ = 0. Variation with respect to Aa yields the Maxwell 
equations with a current of the form 



f = Vpipu'^qf-pu^q'^). 



(25) 



When a vector is found to give the current (25), the Lorenz force term in the MHD-Euler 
equation 

u^{d{hu))p^ = ^F^pf (26) 
P 

can be absorbed into the left hand side, bringing eq. (26) in the canonical form 

u^ojpa = 0, (27) 

where 

■= {dw)a0, Wa := hUa + Va, Va ■= Foc^q'^ , (28) 

because of a relation, 

-pFe^pf = [^q{pU^) + PU^VW] = {dl])apU^ . (29) 

The one-from Wa and its exterior derivative ojaji can be respectively regarded as a canonical 
momentum of a magneto-fluid element and a generalized vorticity. The MHD-Euler equation 
(27) in this case implies conservation of circulation for the magnetized flow, because its curl, 
with do; = 0, yields 

£«^a/3 = 0, (30) 
and the Lagrangian perturbation of this conservation law has a first integral 

Aa;«;3 = 0. (31) 
Substituting the current (25) to the first law (14), we have 

SQ = j !^ ^t^dS + ^~^^^ AdMB + v'^AdCa - v^q^/^Fp^dM^ - {j'^k^ - /k'^)AA0 dSa | 
where we introduced the notation 



(32) 



dMB := pu'^dSa, dS := sdMe, dCa := {hua + r]a)dMB. (33) 

Applying the first law (32) to a sequence of solutions along which the quantities are conserved 
as in Eq. (23) and the circulation of magnetized flow is conserved as in Eq. (31), the first law 
becomes 

SQ = 0, (34) 
or for asymptotically flat systems, such as in the post-Newtonian approximation, dQ = 

SM - nsj = 0. 



3. A formulation for computing magnetized neutron star equilibria 
To compute solutions of single or binary compact objects in equilibrium, we apply a finite 
difference scheme, or a pseudo-spectral method, to a system of basic equations and numerically 
solve the system. Those equations include the Einstein equation, the Maxwell equations, the 
MHD-Euler equation, and the baryon mass conservation equation. Here, we assume that the 
flow is homentropic, and hence the neutron-star matter is described by a one-parameter EOS. 

For the gravitational field, the Isenberg - Wilson - Mathews formulation, the waveless 
formulation, or a full set of Einstein's equations assuming helical symmetry may be used. These 



formulations are based on a 3+1 decomposition of spacetime, and reliable numerical methods 
to solve these equations have been developed [16-18]. We expect that the Maxwell equations 
can be solved using analogous formulations and applying one of the above numerical methods. 

However, for the MHD-Euler equation, when stationarity or helical symmetry is imposed, it 
is no longer an evolution equation, and as a result it is difficult to integrate numerically. In 
self-consistent field methods [16-19], an equilibrium solution is computed using a first integral 
of the (MHD-) Euler equation which exists when the flow is assumed to be either corotational or 
irrotational. Therefore, finding the first integral of the MHD-Euler equation is a key, and also a 
restriction, for computing equilibrium solutions considered in the previous section successfully. 

If we assume that the orbit of the stars is closed and quasi-circular, then the system appears 
stationary in a frame rotating with frequency equal to the orbital frequency fi, and there exists 
an approximate helical Killing vector A;", given by eq. (3), that Lie-derives all variables in (6). 
The Cartan identity implies 

k^{dw)pa = ^-kWa - Vaik'^w-y) (35) 
It should be noticed that £kWa / in general, since is not an observable quantity. One has 

£kWa = ^krja = Fafj£kq'^. (36) 

(In fact, if one assumes £kWa = for a corotational flow, no Lorenz force is exerted on the 
matter). As shown below (see also [20]), the Cartan identity (35) leads quickly to conserved 
quantities when certain assumptions for the flow are made. 



3.1. Irrotational magneto-flow 

If one assumes that the magneto-flow is irrotational, i.e. described by a potential 

Wa = hua + r]a = V^^, (37) 

then the vorticity oj = dw vanishes and the MHD-Euler equation (27) is automatically satisfied 
while the left hand side of eq. (35) vanishes as well, so that 

£fef«a - Vaik'^Wj) = (38) 

The above equation has a first integral iff the term £kWa is the gradient of some scalar function. 
Because of eq. (36), this integrability condition is written as £kql^Fpa = ^af, or, using the 
Cartan identity, 

£-[k,q]^aP = 0) (39) 

where [k,q]" := £kq'^. Then, eq. (35) implies that the quantity 

k'^w^ + f = £ 

is constant throughout the magneto-fluid. 

In this case, we have four variables for the matter, two thermodynamic variables and two 
variables for the velocity fields, {h,p, u*}. The set of equations to solve for these four variables 
is supplemented by a one-parameter EOS, p = p{h), and the normalization of the 4- velocity 
Uau" = —1. The velocity potential $ is obtained from the rest mass conservation equation 
which is an elliptic equation for $ with Neumann boundary conditions on the stellar surface. 
An additional degree of freedom may be obtained from the first integral (40), if it exists. 

It is not trivial to find a vector (7" that satisfies both the integrability condition (39) and the 
ideal MHD condition (20); is not freely specifiable. Instead of solving for a that satisfies 
both conditions, let us assume that 

= (41) 



(40) 



with Zj^q* regarded a function of Aak°'. If we select q"' as above, at least instantaneously on an 
initial hypersurface E, then the first integral is written 

k'^Wa + J Ekq'diA^k'') = S (42) 

It is likely that solutions with neutron stars and strong magnetic fields computed from Eq. (42) 
may not be strictly in equilibrium. Nonetheless, such solutions with strong magnetic fields 
that satisfy a set of hydrostationary equations will serve as interesting initial data sets for 
numerical relativity merger simulations and for further studies of such strongly magnetized 
compact binaries. 



3.2. Corotational magneto-flows 

If the flow is corotational, as is the case for a tidally locked binary or a rigidly rotating non- 
axisymmetric magnetar, then the velocity field can be written it terms of the helical Killing 
vector (3) as 

= (43) 

In this case, the MHD-Euler equation (27) implies that the left-hand side of the Cartan identity 
(35) vanishes, so that eq. (38) again holds. The integrability condition of the latter equation 
then is formally identical to eq. (39) and the first integral of eq. (38) again has the form of eq. 
(40): 

k-tw^ + f = £ (44) 
In the corotational case, the perfect MHD condition (20) also admits a first integral: 

k^A^ = 

This follows from the Cartan identity, 

kl^{dA)0^ = £kAa - Va(A;^A), (46) 

along with the facts that the term kl^{dA)0a vanishes by virtue of eqs. (43) and (20), and that 
the term £kAa can be made to vanish by a gauge transformation in A^ (c.f. Appendix C in 
[22]). One can then incorporate the first integrals (44) and (45) in a self-consistent iteration 
scheme analogous to that outlined for the irrotational case above, coupled to the Einstein- 
Maxwell equations, to construct equilibrium models of rigidly rotating configurations. A scheme 
based on cq. (44), with / = 0, has been implemented in [23] for constructing non-magnetized 
triaxial (Jacobi) spheroids in general relativity and can be extended to model magnetars that 
are triaxially deformed by their magnetic field. 



(45) 



4. Discussion 

The thermodynamic laws derived for helically symmetric Einstein-Maxwell spacetimes with 
magnetized matter can be applied to a non-axisymmetirc single rotating star, or an axisymmetric 
rotating star beyond the stationary, axisymmetric and circular spacetime. Our project involves 
computing quasi-equilibrium configurations of magnetized compact objects with numerical codes 
based on self-consistent field methods. Such methods can be used to model not only binary 
compact objects, but also magnetars, as well as proto neutron stars which are likely to have 
strong magnetic fields. 
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